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ABSTRACT 

We present simulations modeling closed regions of the solar corona threaded by a strong magnetic field 
where localized photospheric vortical motions twist the coronal field lines. The linear and nonlinear 
dynamics are investigated in the reduced magnetohydrodynamic regime in Cartesian geometry. Ini- 
tially the magnetic field lines get twisted and the system becomes unstable to the internal kink mode, 
confirming and extending previous results. As typical in this kind of investigations, where initial 
conditions implement smooth fields and flux- tubes, we have neglected fluctuations and the fields are 
laminar until the instability sets in. But previous investigations indicate that fluctuations, excited by 
photospheric motions and coronal dynamics, are naturally present at all scales in the coronal fields. 
Thus, in order to understand the effect of a photospheric vortex on a more realistic corona, we continue 
the simulations after kink instability sets in, when turbulent fluctuations have already developed in 
the corona. In the nonlinear stage the system never returns to the simple initial state with ordered 
twisted field lines, and kink instability does not occur again. Nevertheless field lines get twisted, but in 
a disordered way, and energy accumulates at large scales through an inverse cascade. This energy can 
subsequently be released in micro-flares or larger flares, when interaction with neighboring structures 
occurs or via other mechanisms. The impact on coronal dynamics and CMEs initiation is discussed. 
Keywords: magnetohydrodynamics (MHD) — Sun: corona — Sun: coronal mass ejections (CMEs) 
— Sun: magnetic topology — turbulence 
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1. INTRODUCTION 

Photospheric convection and the coronal magnetic field 
play a key role in heating the solar corona. For the mag- 
netically confined (closed) regions of the corona it has 
long been suggested that small heating events, dubbed 
nanoflares, continuously deposit energy at the small 
scales and ca n contribute a substanti al fraction of the 
total heating (jParkerll 19721 [l988l Il99l . 

The slow photospheric motions (with a timescale r c ~ 8 
minutes, magnitude v c ~ 1 /cms -1 , and correlation scale 
i c r^j 10 3 km) transfer energy from the photosphere into 
the corona shuffling the footpoints of the magnetic field 
lines. The work done by the denser photospheric plasma 
on the magnetic field lines footpoints injects energy into 
the corona, mostly as magnetic energy. 

The perturbations generated at the photospheric level 
propagate along the loop at the Alfven speed. In a coro- 
nal loop the Alfven velocity associated to the strong 
axial magnetic field Bo is va — Bo/^Airpo ~ 2 x 
10 3 kms~ 1 (po is the average mass density), and con- 
sidering a typical loop length L ~ 4 x 10 4 km we ob- 
tain for the Alfven crossing time ta = L/va ~ 20 s. 
The crossing time is therefore significantly smaller than 
the photospheric timescale: ta « r c . Furthermore as 
magnetic pressure largely exceeds plasma pressure the 
plasma /3 is small (<< 1). Because of the fast Alfven 
timescale and low beta the dynamics of the magneti- 
cally confined solar corona are generally approximated 
as a quasi-static evolution through a sequence of equilib- 
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no, with instabilities leading the system fro m an equilib- 
rium to the next one through r elaxation (Ta vlorl Tl974 
119861 : iHevvaerts fc Priestl 119841 ). Within this frame- 
work many works have studied the relaxation dy nam- 
ics in detail (e.g., lYeates et all 120101 ; iPontin et"aLll2011b 
iWilmot-Smith et al.l[20TTI) . 

On the other hand it has been shown that in some 
instances the validity of that picture is not valid, e.g., 
in reduced magnetohydrodynamic (MH D) simulations of 
the Parker model for coronal heating ( Rappazzo et all 
[20071 [20081 : iRappazzo. Velli fc Einaudll boiOl ). In fact 
that approximation is attained neglecting the velocity 
and plasma pressure in the MHD equations whose solu- 
tion is then bound to be a static force-free equilibrium. 
But the self-consistent evolution of the plasma pressure 
and velocity, although small compared with the domi- 
nant axial magnetic field Bq (and therefore very close in 
first approximation to the equilibrium solution of a uni- 
form axial magnetic field), does not bind the system to 
force- free equilibria and allows the development of turbu- 
lent nonlinear dynamics with formation of field-aligned 
current sheets where a significant heating occurs. 

Most numerical simulations of a simple model coro- 
nal loop in Cartesian geometry threaded by a strong 
magnetic field shuffled at the top and bottom plates by 
photospheric motions have used as boundary velocity 
an incompressi ble field with all wavenumbers of order 
~ 4 ex cited (|Einaudi et al.l 119961 : I Dmitruk fc Gomez l 
19971: [Georgoulis et al.l 119981; IDmitruk fc Gomezl 
1999; Einaudi & Velli 1999; Dmitruk et al. 2003; 



Rappazzo et al. 2007, 2008; Rappazzo & Velli 2011): in 
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real space this corresponds to distorted vortices with 
length-scale ~ 10 3 km one next to the other filling 
the whole photospheric plane (Rap pazzo et al.l |.2008, 
Figures 1 and 2). This configuration does not give rise to 
instabilities. The system transitions smoothly from the 
linear to the nonlinear stage where integrated physical 
quantities like energies and dissipation fluctuate around 
a mean value in what is best described as a statistically 
steady state and the energy deposited at the small 
scales is approximately in the nanoflare range for the 
numerous heating events. 

This disordered vortical forcing mimics a uniform and 
homogeneous convection and the resulting coronal dy- 
namics give rise to a basal background coronal heating 
within the lower limit of the observational constraint. 

With this kind of forcing the system is not able to ac- 
cumulate a significant amount of magnetic energy to be 
subsequently released in more substantial heating events 
like microflares and flares. It is therefore pivotal to im- 
plement different kinds of photospheric forcings to under- 
stand the role of convective motions in coronal heating 
and the physical mechanisms and conditions for a signif- 
icant storage of magnetic energy and its release. 

Shearing or twisting the field lines might appear 
as the most straightforward way to make the sys- 
tem accumulate energy. To thi s end, in recent work 
(Rappaz zo. Velli fc Einaudlll2010f ) we have implemented 
a ID velocity forcing with a sinusoidal shear flow at 
the boundary, with v y (x) oc sin(27rkx/£), wavenumber 
k = 4, spanning the whole photospheric plane {I is the 
cross-section length). Initially the magnetic field that 
develops in the coronal loop is a simple map of the pho- 
tospheric velocity field, i.e., b y (x) oc (t/TA)sin(2irkx/£), 
with its intensity growing linearly in time. A sheared 
mag netic field is know n to be subject to tearing instabil- 
ity j^rtF^F^l][T963), in fact magnetic energy accumu- 
lates until a tearing instability sets in, magnetic energy 
is released and the system transitions to the nonlinear 
stage. 

On the other hand continuing the simulation we found 
that, once the system has become fully nonlinear the dy- 
namics are fundamentally different: the magnetic shear is 
not recreated. Once fluctuations are present, the orthog- 
onal magnetic field (b x and b y ) is organized in magnetic 
islands with X and O points, the nonlinear terms do not 
vanish any more and energy can be transported efficiently 
from large to small scales where it is dissipated. 

In the nonlinear stage the dynamics are very similar 
whether the forcing velocity is a shear flow or made of 
disordered vortices, and magnetic energy is not stored 
efficiently, therefore larger releases of energy are not pos- 
sible. A shear photospheric flow can give rise to a sheared 
coronal magnetic field only in the unlikely condition that 
no relevant perturbations are present in the corona, or for 
very strong shear flows with velocity higher than the typ- 
ical photospheric velocity of 1 /cms -1 , when the pertur- 
bations naturally present in the corona can be neglected. 

All photospheric forcings described so far fill the entire 
photospheric plane. Spatially localized velocity fields, 
like a single vortex that does not fill the entire plane, 
might be able to induce a higher storage of magnetic 
energy in the corona. 

A vortex twists the magnetic field lines, and the result- 
ing helical magnetic structure is kink unstable and widely 



used to model co ronal loops (Batv & Hevvaerts 1996; 
Velli et all 119971: iLionello et al.l 119981 : iBrowning et~aT1 
20081 : iHood eTal]l2QQ9l ). In these studies a twisted mag- 
netic field lines structure is used as initial condition, but 
this is assumed to have been induced by photospheric 
motions shuffling their footpoints, while the field lines 
are actually line-tied to a motionless photosphere. 

To understand the dynamics of this photospherically 
driven system we have performed numerical simulations 
appl ying localized vor ti ces at the photosp heric bound- 
ary. IMikic et all (|1990D : iGerrard et al.l (|2002f ) performed 
boundary forced simulations, but they stop just after 
kink instability sets in. In this paper we continue the sim- 
ulations for longer times. This allows us to understand 
the dynamics of the system both when initially only an 
axial uniform magnetic field is present and a smooth or- 
dered flux-tube with twisted field lines gets formed and 
kink instability sets in, and the later dynamics when the 
localized boundary vortex twists the footpoints of a dis- 
ordered magnetic field where magnetic fluctuations and 
small scales are already present and field lines are no 
longer smooth. 

Furthermore, highly twisted magnetic structures, such 
as flux ropes, are broadly used to model the i nitia- 
tion o f coronal mass eject io ns (CMEs, e. g.. see lLow 
(2001b; lAmari fc Alyl (|2009h : IChenl (|201lh : frorok et al. 
(j2011f h and references therein). To advance our un- 
derstanding of eruption initiation it is therefore impor- 
tant to understand if and under which conditions photo- 
spheric motions can self-consistently generate flux ropes 
(jAmari et al.lll999l ). or if these structures can only be 
advected into the corona from sub-photospheric regions 
via emerging flux. 

The paper is organized as follows. In § [2] we describe 
the basic governing equations and boundary conditions, 
as well as the numerical code used to integrate them. In 
§ [3] we discuss the initial conditions for our simulations 
and briefly summarize the linear stage dyna mics more 
extensively detailed in Rapp azzo et al.l ([20 08) . The re- 
sults of our numerical simulations are presented in § HJ 
while the final section is devoted to our conclusions and 
discussion of the impact of this work on coronal physics. 



2. GOVERNING EQUATIONS 

We model a coronal loop as an axially elongated Carte- 
sian box with an orthogonal cross section of size £ and 
an axial length L embedded in an homogeneous and uni- 
form axial magnetic field Bq = Boe z aligned along the 
^-direction. Any curvature effect is neglected. 

The dyn amics are integrated with the equations 
of RMHD (iKadomtsev fc Pogutsd 119741 : iStraussl fl976L 
lMontgomervlll982f ). well suited for a plasma embedded 
in a strong axial magnetic field. In dimensionless form 
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they are given by: 

d\x | , ^ x ^ ( b . 2 

- g A + („ i .VJu i =-vJ P +-±- 



(b x • V x )b x +ca 



-gf + K • V,) b ± = (b ± -V ± )u ± +c A -^ 

V x • u x = 0, V x • b x = 0, (3) 

where u x and b x are the velocity and magnetic fields 
components orthogonal to the axial field and p is the 
plasma pressure. The gradient operator has components 
only in the perpendicular x-y planes 



± ~ ex dx~ + ev dy~ 



(4) 



while the linear term oc d z couples the planes along 
the axial direction through a wave-like propagation at 
the Alfven speed ca- Incompressibility in RMHD equa- 
tions f ollows from t he large value of the axial magnetic 
fields (lStrausslll976|) and they remain valid also for low 



P sys tems ([Zank fc Matthaeus 1992; Bhattachariee et al 
1998) such as the corona. Furthermore iDahlburg et al 
(|2012[ ) have recently performed fully compressible simu- 
lations of a similar Cartesian coronal loop model, show- 
ing that the inclusion of thermal conductivity and radia- 
tive losses, that transport away the heat produced by the 
small scale energy dissipation, keep the dynamics in the 
RMHD regime. 

To render the equations nondimensional, we have first 
expressed the magnetic field as an Alfven velocity [b — > 
bj y / 47rpo], where po is the density supposed homogeneous 
and constant, and then all velocities have been normal- 
ized to the velocity u* = 1 /cms -1 , the order of magni- 
tude of photospheric convective motions. 

Lengths and times are expressed in units of the per- 
pendicular length of the computational box £* = i and 
its related crossing time t* = £*/u*. As a result, the 
linear terms ex d z are multiplied by the dimensionless 
Alfven velocity ca = va/u*, where va = Bo/y/Airpo is 
the Alfven velocity associated with the axial magnetic 
field. We use a computational box with an aspect ratio 
of 10, which spans 



0<x,y<l, 



< z < 10. 



(5) 



Our forcing velocities have a linear scale of ~ 1/4 that 
corresponds to the convective scale of ~ 1,000 km in 
conventional units, thus the box extends (4,000 km) 2 x 
40, 000 /cm. 

The index n in the diffusive terms ([I])-© is called dis- 
sipativity and for n > 1 these co rrespond to so-called 
hyperdiffusion (e.g., Biskamp 2003). For n = 1 standard 
diffusion (Re 1 = Re) is recovered and in this case the 
kinetic and magnetic Reynolds numbers are given by: 

Re = , Re m = r — , (6) 

V T]C Z 

where c is the speed of light, and numerically they are 
given the same value Re = Re m . 



Table 1 

Simulations summary 



Run 


n x 


x n y 


X 


n z 


Re 4 


t"max 


A 


512 


x 512 


X 


208 


1 x 10 19 


1,900t a 


B 


256 


x 256 


X 


192 


5 x 10 16 


ll,000r A 



Note. — The numerical grid resolution is n x x n y x n z . The 
next columns indicate respectively the value of the hyperdiffusion 
coefficient Re 4 and the simulation time span. 



In the simulations presented in this paper we use hy- 
perdiffusion with n = 4. Hyperdiffusion is used be- 
cause the implemented boundary velocity forcings and 
the magnetic flux tubes induced initially are localized to 
a small area of the computational box, and the dynamics 
would be dramatically diffusive with standard diffusion 
at a reasonable resolution (see next section §[3] for a more 
detailed discussion). 

Our parallel code RMH3D solves numerically Equa- 
tions (PO)-© written in terms of the potentials of the 
orthogonal velocity and magnetic fields in Fourier space, 
i.e., we advance the Fourier components in the x- and y- 
directions of the scalar potentials. Along the z-direction, 
no Fourier transform is performed so that we can im- 
pose non-periodic boundary conditions (§[3|), and a cen- 
tral second-order finite-difference scheme is used. In the 
x-y plane, a Fourier pseudospectral method is imple- 
mented. Time is discretized with a third-order Runge- 
Kutta method. F or a more detailed des cription of the 
numerical code see iRappazzo et ahl (j2007l [2008). 

3. BOUNDARY CONDITIONS AND LINEAR 
STAGE DYNAMICS 

Magnetic field lines are line-tied to the top and bottom 
plates (z = and L) that represent the photospheric sur- 
faces. Here we impose, as boundary condition, a velocity 
field that convects the footpoints of the magnetic field 
lines. Along the x and y directions periodic boundary 
conditions are implemented. 

All simulations (see Table [1]) employ a circular vor- 
tex applied at the top plate z = L. The velocity po- 
tential for this vortex is centered in the interval / = 
[1/2 - 1/8, 1/2 + 1/8] of linear extent 4 = 1/4 and van- 
ishes outside: 



ip(x,y) — ■= sin 2 

for x,y G /, and 



47T ( x — — 



(p = for x,y £ I. 



The velocity is linked to the potential by \x ± 
and its components are: 



(7) 

Vip x e z 



L / \ ,2. 2 

u x (x,y) = +-j= sm 

L / \ 2 . 

u y (x,y) = --^sm 



47r ( x — i 



8tt [x - - 



87r [y - o 



4tt [y - - 



(8) 
(9) 



in the interval / and vanish outside. As shown in Fig- 
ure [1] Equations jl])-(|9j) describe a counter-clockwise vor- 
tex centered in the middle of the plane z = L and 
has circular streamlines, with a slight departure from 
a perfectly circular shape toward the edge of the inter- 
val /. Averaging over the surface / the velocity rms is 
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Figure 1. Circular vortex employed as boundary velocity forcing 
in the presented simulations. Above: streamlines and profile of its 
absolute value \u\. Below: plot of the velocity y-component as a 
function of x at y = 0.5. 
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fields used in our previous works ([R appazzo 



1/2, the same value of the boundary velocity 

et al.jj2 007. 



apt 

2008: lRappazzo. Velli fc Einaud i 2010; Rappaz zofc Vellil 
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In all simulations a vanishing velocity is imposed at 
the bottom plate z = 0: 



u° (x,y)=0. 



(10) 



At time t = we start our simulations with a uniform 
and homogeneous magnetic field along the axial direction 
Bo = Boe z . The orthogonal component of the velocity 
and magnetic fields are zero inside our computational 
box u_l = b_L =0, while at the top and bottom planes 
the vortical velocity forcing is applied and kept constant 
in time. 

We briefly summarize and specialize to the case con- 
sidered in this p aper the linear stage a nalysis described 
in more detail in Rappaz zo et al.l ([20081 ). In general for 
an initial interval of time smaller than the nonlinear 
timescale t < r n i, nonlinear terms in Equations (pp)- ((3]) 
can be neglected and the equations linearized. For sim- 
plicity we will at first neglect also the diffusive terms 
and consider their effect later in this section. The so- 
lution during the linear stage with a generic boundary 
velocity forcing u L , and u° = 0, (respectively at the top 
and bottom planes z = L and 0) is given by: 



b ± (x,y,z,t) = — u L , 

TA 

u ± (x,y,z,t) = -u L . 



(11) 
(12) 



where ta = L/va is the Alfven crossing time along the 
axial direction z. The magnetic field grows linearly in 
time, while the velocity field is stationary and the order 
of magnitude of its rms is determined by the boundary 
velocity profile. Both are a mapping of the boundary 
velocity field u L . 

For a generic forcing the solution (fTT]) - ([T2]) is valid only 
during the linear stage, while for t > r n \ when the fields 
are big enough the nonlinear terms cannot be neglected. 
Nevertheless there is a singular subset of velocity forcing 
patterns for which the generated coronal fields ([TT]) -([T2 ]) 



have a vanishing Lorentz force and the nonlinear terms 
vanish exactly. 

This subset of patterns is characteriz e by having the 
vorti city constant along the streamlines ([Rappazzo et al.l 
2008). In this case magnetic energy grows quadratically 
in time until some instability eventually sets in. Two 
kind of velocity patterns can be identified: a) ID pat- 
terns with their streamlines all parallel to each other, 
like a shear flow, or b) a radial pattern with circular 
streamlines, like a circular vortex. 

Since in the linear stage the coronal fields are a map- 
ping of the boundary velocity (fTT]) - ([T2]) . a shear flow in- 
duces ash^aredm^^ tearing instabil- 
ities (Rappa zzo. Velli fc Einau di 2010), while the vorti- 
cal flows considered in this paper twist the field lines into 
helices subject to kink instabilities. The vortex (j8j)-(j9j) 
is not perfectly circular as the streamlines depart from 
an exact roun d sh ape toward the edge (Figure HJ, but as 
we show in § 14.11 field line tension adjusts the induced 
coronal orthogonal field lines in a round shape. 

So far we have neglected the diffusive terms in the 
RMHD Equations (P)-(j3j). In § H we show that for this 
kind of problem the use of hyperdiffusion is crucial, oth- 
erwise the dynamics are dominated by diffusion. Over- 
looking t his numerical fact can resu lt in misleading con- 
clusions (jKlimchuk et a l. 2009, 2010), upon which we will 
comment in § |4j 

Here we want to understand the diffusive effects on the 
linear dynamics, i.e., when nonlinear terms are negligible 
or artificially suppressed by a low numerical resolution. 
We now consider the effect of standard diffusion (case 
n = 1 in eqs. ([I])-®) on the solutions ([TT]) - ([T2]) : these are 
the solutions of the linearized equations obtained from 
(PQ)-(j2]) retaining also the diffusive terms. 

In the linear regime, as the magnetic field grows in 
time ([TT]h the diffusive term (V^ h ± ocb ± /£ 2 ) becomes 
increasingly bigger until diffusion balances the magnetic 
field growth, and the system reaches a saturated equilib- 
rium state. Including diffusion the magnetic field evolves 
as 



b ± (x,y,z,t) = u (x,y) — 

TA 



exp ( 

TR 



, (13) 



i.e., for times smaller than the diffusive timescale tr 
Equation ([TT]) is recovered with the field growing linearly 
in time, while for times bigger than tr the field asymp- 
totes to its saturation value. The diffusive timescale asso- 
ciated with the Reynolds number Re is tr = £ 2 Re/ (2tt) 2 
where £ c is the length-scale of the forcing pattern, that 
for the pattern (O-® is given by t c ~ £/4 where t is the 
orthogonal computational box length. 

The total magnetic energy Em and ohmic dissipation 
rate J will then be given by 



E 



1 

M = - 



dx b 2 



J M 



1 I 1 

1 — exp 

TR 



J - — f dx \ 2 - J sat 
J ~ Re J V &X} ~ J 



1 — exp 

TR 



, (14) 
(15) 



For times smaller than the diffusive timescale tr both 
quantities grow quadratically in time, while for t > 2tr 
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Figure 2. Run A: Magnetic (Em) and kinetic (Ek) energies as 
a function of time (ta = L/va is the axial Alfven crossing time). 
The dashed curves show the time evo lution of magnetic energy if 
the system were unperturbed [ea.J14))]. or nonlinearity suppressed, 
e.g., by numerical diffusion (see The case R = oo corresponds 
to the linear case with no diffusion, attained (at the large scales) 
with the implementation of hyperdiffusion. 

they asymptote to their saturation value Effi and J sat - 



rrisat 



i%c\Re 2 
2L(2n) 4 



J s 



O ipsat 
f M 

(4ir) 2 £ 2 c Re' 



(16) 



Magnetic energy saturates to a value proportional to the 
square of both the Reynolds number and the Alfven ve- 
locity, while the heating rate saturates to a value that is 
proportional to the Reynolds number and the square of 
the axial Alfven velocity. 

Even though we use grids with ~ 512 2 points in the 
x-y plane, the timescales associated with ordinary diffu- 
sion are small enough to affect the large-scale dynamics, 
inhibiting the development of instabilities and nonlinear- 
ity. The diffusive time r n at the scale A associated with 
the dissipative terms used in Equations ([I])-© is given 

by 

r n ~ Re n A 2n . (17) 

For n = 1 the diffusive time decreases relatively slowly 
toward smaller scales, while for n = 4 it decreases far 
more rapidly. As a result for n = 4 we have longer 
diffusive timescales at large spatial scales and diffusive 
timescales similar to the case with n = 1 at the resolu- 
tion scale. Numerically we require the diffusion time at 
the resolution scale Xmin — where N is the number 
of grid points, to be of the same order of magnitude for 
both normal and hyper-diffusion, i.e., 



Re 1 



7V 2n 



Re 



Re, TV 



2(n-l) 



(18) 



Then for a numerical grid with N = 512 points that 
requires a Reynolds number Re 1 = 800 with ordinary 
diffusion we can implement Re 4 ~ 10 19 (table [TJ, remov- 
ing diffusive effects at the large scales and allowing, if 
present, the development of kink instabilities and non- 
linear dynamics. 

4. NUMERICAL SIMULATIONS 

In this section we present the results of the numerical 
simulations summarized in Table [U Simulations A and 
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Figure 3. Run A: Ohmic (J) dissipation rate and the integrated 
Poynting flux S (the injected power) versus time. Viscous dissi- 
pation is negligible respect to the ohmic contribution. Inset shows 
the ohmic dissipative peak corresponding to the development of 
kink instability. 

B have the same parameters, but simulation B employs 
a lower resolution to achieve a very long duration. In 
all simulations the vortical velocity pattern ([I])-(j9j) is ap- 
plied at the top plate z = L, and a vanishing velocity at 
the bottom plate z = 0. Initially no perpendicular mag- 
netic or velocity field is present inside the computational 
box = u ± =0, and the system is threaded only by 
the constant and uniform field Bo = B$e z . The compu- 
tational box has an aspect ratio of 10, with £ = 1 and 
L = 10. 

4.1. Run A 

We present here the results of run A, a simulation per- 
formed with a numerical grid of 512 x 512 x 208 points, 
and hyperdiffusion coefficient Re^ = 10 19 with diffusivity 
n = 4. The Alfven velocity is va = 200kms _1 , corre- 
sponding to a nondimensional ratio ca — va/u* = 200. 
The total duration is ~ 1, 900 axial Alfven crossing times 
t a = L/v A . 

Figures [2]|3] show the temporal evolution of the total 
magnetic and kinetic energies 



E 



A I 



\ idVb 2 



E 



K 



\ IdVv? 



the total ohmic dissipation rate 



J Re 



(19) 



(20) 



and 5, the power injected from the boundary by the work 
done by convective motions on the field lines' footpoints 
(see Equation (|22]) ), along with some saturation curves 
for magnetic energy (fT4|) . Additionally Figure H] shows 
snapshots of the magnetic field lines of the orthogonal 
component and electric current j = j z , the lea ding 
order component in RMHD ordering (Strauss 1976), at 
selected times in the mid-plane z = 5. 

The circular vortical velocity field dH])-© applied at 
the top boundary (z = 10) initially induces velocity and 
magnetic fields in the computational box that follow the 
linear behavior given by Equations (fTT]) - (fT2]) . i.e., they 
are a mapping of the velocity at the boundary with the 
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magnetic field increasing linearly in time (Figure 2J times 
t = 0.61 ta and 80.64 ta). In the linear stage (t < 83 ta) 
magnetic energy is well-fitted (Figure [2j) by the linear 
curve ([T4]) in the limit Re — >• oo, i.e., in the absence of 
diffusion (indeed in this limit the curve can be obtained 
directly from the linear Equation ([TT]) ). This is because 
we are using hyperdiffusion that effectively gets rid of 
diffusion at the large scales. 

Two other magnetic energy linear diffusive saturation 
curves are drawn for Re = 800 and 400, typical Reynolds 
numbers used in our previous simulations with standard 
diffusion n = 1 and orthogonal grids with respectively 
512 2 and 256 2 grid points (see, e.g., Rapp azzo et al.l 
2008). Their saturation level is very low compared to 
the magnetic energy values when kink instability devel- 
ops (t ~ 83 ta) and in the following nonlinear stage. This 
is because the vortex and induced magnetic field occupy 
only a limited volume elongated along z at the center of 
the x-y plane: at these scales diffusion dominates with 
these resolutions using standard diffusion. 

For this reason the use of hyperdiffusion is crucial to 
study this problem, otherwise diffusion dominates and a 
balance between the injection of energy from the bound- 
ary and its numerical removal by diffusion is reached 
very soon, inhibiting the development of kink instabil- 
ity and nonlinear dynamics. This diffusi ve linear regime 
was reached in previous simulations by iKlimchuk et al.l 
(20091. l201(I ). where four similar vortices were applied at 
the boundary. Therefore their conclusion that nonlin- 
ear dynamics or instabilities (not to mention turbulence) 
cannot develop in such physical systems is simply a nu- 
merical issue: this can be overcome adopting hyperdiffu- 
sion as we have done here or, alternatively, implementing 
grids with much higher resolutions that require imprac- 
tically large numerical resources. 

The localized boundary vortex (shown with a colored 
contour in Figure [5} generates a mostly poloidal mag- 
netic field confined to the axial volume in correspondence 
of the vortex, resulting in helical field lines for the total 
magnetic field (Figure [5j time t = 60.57 ta)- Outside 
this volume the poloidal field vanishes and only the axial 
field Bo is present. Ampere's law then guarantees that 
the total net current is zero. As shown in Figured] in the 
linear stage (t = 0.61 ta, 80.64 ta) there is a stronger up- 
flowing current concentrated in the middle, and a weaker 
ring-shaped down-flowing current distributed at the edge 
of the flux-tube. 

This magnetic configuration is well known to be kink 
unstable, and is simila r to the NC (Null Cu rrent) force- 
free model studied bv iLionello et al.l (jl998[ ). The main 
differences are that their axial field Bo is not uniform, 
dropping by ~ 50% outside the flux tube, and that the 
field lines are line-tied to a motionless photosphere. They 
performed a linear stability analysis of this configuration 
finding that there is a critical axial loop length L cr u be- 
yond which the system is unstable and has a constant 
growth rate jta ~ 0.02. They also examined other equi- 
libria with net current finding a similar qualitative be- 
havior, with variations for the critical length and growth 
r ates. 

ILionello et "all (|1998D found that for the NC case the 
ratio of the axial critical length over the cross-length of 
the flux-tube is L cr i t /£ c ~ 9. In the case considered here 
the ratio of the axial length (L = 10) over the cross- 



length of the flux tube (the extent of the boundary vor- 
tex £ c = 1/4, Equation (O) is L/£ c = 40, therefore it is 
fully in the unstable region. Of course at a given length 
(beyond the critical length) there is also a critical twist 
beyond which the configuration is unstable. In our sim- 
ulations the system is continuously forced at the bound- 
ary, and in the linear stage the twist grows linearly in 
time (from Equation ([TT]) . as the twist is proportional to 
b±/Bo), thus such a critical twist is certainly attained. 

In our case the "equilibrium" solution is not static but 
is given by the linear solution ([TT]) . indicated here with 
bii n , with the magnetic field growing linearly in time 
while mapping the boundary vortex. Thus we compute 
the perturbed magnetic energy as 

E M = [ d% |b-b, m | 2 . (21) 
Jv 

We find that in the linear stage this quantity grows ex- 
ponentially in time, obtaining for the p erturbed mag- 
netic field a growth rate 7TA ~ 0.02, as ILionello et al.l 
(1998) for their NC equilibrium model. This growth rate 
is also confirmed by the fact that kink instability sets 
in at t ~ 83 ta (Figures [2] and [3]) and I/7 ~ 50 ta- As 
mentioned in § [3] the forcing boundary vortex departs 
from an exact circular shape at its edges where its vor- 
ticity is not exactly constant along the streamlines, thus 
there is a small Lorenz force for the resulting magnetic 
field ([TT]) . This small difference in the linear field acts as 
a perturbation. 

Additionally ILionello et al.l (jl998[ ) found out that con- 
figurations with zero net current are unstable to the in- 
ternal kink mode (opposed to the global kink mode for 
configurations with a net current), for which magnetic 
perturbations and the radial displacement of the plasma 
column are confined within the original flux tube. This 
is found also in our simulation as shown in Figure |4] at 
the onset of the nonlinear stage at t = 83.85 ta, when 
the plasma displaces inside the flux tube toward its edge 
where a strong current sheet forms. 

The internal kink mode releases almost 90% of the ac- 
cumulated energy around time t ~ 83.5 ta (Figure [2]) 
in correspondence of the big ohmic dissipative peak 
shown in Figure [3l The released energy is AE ~ 
10 3 x 10 22 erg = 10 25 erg, in the micro-flare range (the 
factor to convert energy into dimensional units, given 
our normalization choice discussed in § [3J is 10 22 , i.e., 
1 — »> 10 22 erg). As a result of the kink instability 
magnetic reconnection occurs (Figure (H t = 85.05 ta) 
and the magnetic field lines get substantially unwind as 
shown in Figured (times t = 60.57 ta and 100.78 ta) 
with field lines twisting only ~ 180° after the instability. 

In summary, during the linear stage, the transition to 
an d the first pha s e of th e nonlinear regime, the analysis 
of ILionello et al.l (jl998[ ) is fully confirmed also for the 
photospherically driven case considered here: the system 
forced by a circular vortex is unstable to an internal kink 
mode, releases most of the stored magnetic energy and 
magnetic rec onnection untwists the field lines. Linear 
calculations (Baty 2001) show that similar dynamics are 
expected also for different configurations with different 
aspect ratios and magnetic guide field values, except for 
those that fall below the instability threshold. 

The phenomenology described so far is also in agree- 
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Figure 4. Run A: Axial component of the current j (in color) and field lines of the orthogonal magnetic field in the midplane (z = 5) at 
selected times covering the linear and nonlinear regimes up to t ~ 600 ta- At the begin ning of the linear stage (t = 0.61 ta) the orthogonal 
magnetic field is a mapping of the boundary vortex [see linear analysis, Equation Still in the linear stage but at later times 

it = 80.64 r^) the field line tension straightens out in a circular shape the vortex mapping. An internal kink mode develops (t ~ 83.85 ta) 
and the instability transitions the system to the nonlinear stage. In the fully nonlinear stage the field lines are still circular, but in a 
disordered way, exhibit a broad range of scales, including current sheets, and steadily occupy a larger fraction of the computational box. 
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Figure 5. Run A: Lateral and top views of magnetic field lines at selected times. In the linear stage (t = 60.57 ta) the boundary vortex 
(shown in color in the plane z = 10) twists into an helix the magnetic field lines in the corresponding region underneath the vortex. Those 
outside this region remain straight, a sample of which is shown in red. Kink instability releases magnetic energy and untwists the field lines 
(t = 100.78 ta), that in the nonlinear stage maintain an approximately constant twist ~ 180°. But with time the region where field lines 
are twisted increases its volume until it fills the whole computational box (t ~ 1211.78 ta)- The box has been rescaled for an improved 
visualization, the axial length (along z) is ten times the length of the orthogonal cross section (along x-y). A lateral view is shown only in 
the first panel, at later times a top view is preferred for a better visualization, as from the side the field lines appear overlapped to each 
other. 
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Figure 6. Run A: Lateral {left column) and top {right column) views of isosurfaces of the squared current j 2 at selected times, respectively 
during the linear stage {t ~ 60.57 r^), right after kink instability (t ~ 100.78 ta), and in the fully nonlinear regime {t ~ 1211.78 ta)- In 
each panel are shown three isosurfaces of j 2 , corresponding respectively to 15% (green), 5% (red) and 2% (blue) of the maximum of j 2 in 
the box at each time. As is typical of current sheets, isosurfaces corresponding to higher values of j 2 are nested inside those corresponding 
to lower values. Although the region where field lines are twisted increases in time (Figure , the current sheets' filling factor remains 
small. The box has been rescaled for an improved visualization as in Figure [5] 
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Figure 7. Run A: Magnetic (Em) and kinetic (Ek) time- 
averaged energy spectra as a function of the orthogonal wavenum- 
ber n ± The inset shows the kinetic energy spectrum of the bound- 
ary velocity vortex [Equations Q-©] applied at the top plate 
z = 10. 



merit with that of three-dimensional simul ations with a 
realistic geometry (jAmari fc Lucianll 12000) . In particu- 
lar strong nonlinearities persist right after the instability 
occurs (t = 85.05 ta and 100.78 ta), when the system 
cannot be described as a constant-a force-free state. An 
inverse cascade of magnetic energy is observed, as the 
orthogonal magnetic field acquires longer scales and the 
overall volume occupied by twisted field lines increases, 
as shown in Figure H] just before (t = 85.05 ta) and af- 
ter (t = 100.78 r^) the instability. In llAmari fc Lucianil 
(2000) this corresponds also to an inverse cascade of mag- 
netic helicity, corresponding in the RMHD case to an in- 
verse cascade of the square potential i/j (see the end of 
this section and our discussion in § 5 for more about this 
quasi- invariant analogous to helicity in RMHD). 

On the other hand, at later times the dynamics are 
certainly surprising when, in the fully nonlinear stage, 
fluctuations created by the kink instability are present 
in the corona. For t > 100 ta magnetic energy increases 
steadily, while kinetic energy remains small (Figure [2]). 
This is in contrast to all our previous simulations with 
space-filling boundary motions, either distorted vort ices 
([Rappazzo et al.ll2007l 120081: [R appazzo & Vellil l201lD or 
shear flows ([Rappazzo. Velli fc Einaudil l2010f ). when in 
the nonlinear regime a magnetically dominated statisti- 
cally steady state was reached where integrated quantities 
would fluctuate around an average value (with velocity 
fluctuations smaller than magnetic fluctuations). 

In our case ohmic dissipation J and the integrated 
Poynting flux S do reach a statistically steady state (Fig- 
ure [3]) . The integrated Pointing flux 

S = c A J dab ± • u L , (22) 

z=L 

is the power entering the system at the boundaries as a 
result of the work done by photospheric motions on the 
footpoints of magnetic field lines (u L is the photospheric 
forcing velocity). But in contrast to our previous results, 
here the power does not balance on the average the dissi- 
pation rate, its average is slightly higher resulting in the 
magnetic energy growth shown in Figure O 
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Figure 8. Run A: Magnetic energy modes versus time. While 
modes with wavenumber n ± > 3 fluctuate around a mean value, the 
first two modes increase steadily showing that an inverse cascade 
occurs. 



In physical space the dynamics are surprising in two 
ways. First, after the kink instability, even though we 
continue to stir the field lines' footpoints with the same 
vortex, no further kin k instability develops. Analogousl y 
to the shear flow case (Rappazzo, Velli & Einaudi 2010), 
once the system transitions to the nonlinear stage the 
magnetic fluctuations generated during the instability 
do not have a vanishing Lorentz force. In fact around 
t ~ 100 r^, at the end of the big dissipative event, the 
topology of the orthogonal component of the magnetic 
field is characterized by circular, but distorted, field lines 
(Figured]). Naturally the Lorentz force does not vanish 
now and the vorticity is not constant along the stream- 
lines. Nonlinear terms do not vanish as they do during 
the linear stage for t < 83 ta- When they vanish mag- 
netic energy can be stored, without getting dissipated, 
into an ordered flux-tube with helical field lines (Fig- 
ure t = 60.57 r^), and matching perfectly round or- 
thogonal magnetic field lines (Figure 2J t = 80.64 ta)- 
But now nonlinearity continuously transfers energy from 
large to small scales where it is dissipated. In physi- 
cal space small scales are not uniformly distributed, but 
they are organized in field- aligned current sheets. These, 
once formed during the onset of the nonlinear stage, per- 
sist throughout the subsequent dynamics (as shown in 
Figures 2] and [6]) , with the energy cascade continuously 
feeding them. Second, the photospherical vortical mo- 
tions do not give rise to an orderly helical flux-tube as 
in the linear stage (Figure [5j t = 60.57 ta)- However, 
magnetic field lines get twisted, but in a disordered way 
(Figure [5] and H t > 100.78 ta)- A new phenomenon 
occurs: on longer timescales the magnetic field acquires 
longer spatial scales (Figured]), the volume where field 
lines are twisted increases (Figure [5]), while the current 
exhibits always a small filling factor occupying a small 
fraction of the volume (Figure [6]) . 

To better understand these phenomena we need to in- 
vestigate the energy dynamics in Fourier space. We con- 
sider the spectra in the orthogonal x-y plane integrated 
along the z direction. As they are isotropic in the Fourier 
k x ~k y plane we compute the integrated ID spectra, so 
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Figure 9. RunB: Over long time-scales also the first two magnetic 
energy modes saturate fluctuating around a mean value. The first 
mode dominates, and the amplitude of its fluctuations corresponds 
to releases of energy of ~ 2 x 10 25 erg, in the micro-flare range. 
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Figure 10. RunB: Ohmic (J) dissipation rate and the integrated 
Poynting flux S (the injected power) versus time. Inset shows the 
ohmic dissipative peak corresponding to the development of kink 
instability. 



that for the total magnetic energy Em we obtain: 
l i 
E M = \ Jdz JJdxdybl 


^ N 

= \ [dz Ib| 2 (k, z) = <T, E M (n), (23) 

k n=l 

where n indicates the shell in /c-space with wavenum- 
ber k = (fc, I) G Z 2 included in the range n — 1 < 
(fc 2 + I 2 ) 1 / 2 < n, and N is the maximum wavenumber 
admitted by the numerical grid (corresponding to the 
smallest resolved orthogonal scale). 

The time averaged magnetic and kinetic energy spec- 
tra as a function of wavenumber are shown in Figure [7J 
the inset shows the spectrum of the boundary vortex' 
kinetic energy (see Equation (O). Photospheric mo- 
tions therefore inject energy at wave numbers between 
2 and 7 (see Equation ([22]) ). the system is magneti- 
cally dominated and the power-laws exhibited at higher 
wavenumber s, in the inert ial range, are similar to those 
obtained with previous space-filling boundary forcings 
(IRappazzo et al.l 120081; IRappazzo. Velli fc Einaudil 120101 : 
iRappazzo fc V elli 201 l|), with the spectrum of magnetic 
energy much steeper than that of kinetic energy. 

However the time- average of the low- wave number 
modes hides an interesting dynamics. Figure [8] shows the 
first five magnetic energy modes as a function of time. 
While modes with wavenumbers n > 3 after the kink 
instability fluctuate around a mean value, the first two 
modes n = 1,2 grow steadily with mode n = 1 becom- 
ing prevalent. This shows that an inverse cascade takes 
place. While the direct cascade transfers energy from 
the injection scale toward small scales (current sheets) 
where energy is dissipated, analogously the inverse cas- 
cade transfers energy toward the large scales (modes 1 
and 2) where no dissipative process is at work and conse- 
quently energy accumulates. In physical space this pro- 
cess gives rise to the large scales that the magnetic field 
acquires in the orthogonal direction, shown in Figures |4] 
and [5j discussed previously. In the RMHD system with 



boundary conditions as we apply here there is no strict 
invariant known to follow an inverse cascade, such as 
magnetic helicity in 3D M HD or the square of the vec- 
tor potential in 2D MHD (jBiskampI 120031 : iBergerl Il997t 
iBrandenburg fc Matthaeusll2004D . RMHD resembles the 
2D MHD case in the sense that though the square of 
the vector potential is not conserved, the terms violating 
conservation arise only from the boundaries in the axial 
direction. A dynamical magnetic inverse cascade mecha- 
nism is therefore still active, impeded only by the inputs 
coming from photospheric motions at the boundary, and 
this explains the accumulation of magnetic energy at the 
largest transverse scales. 

4.2. RunB 

The simulation described in the previous section 
(run A) has a duration of ~ 1200 ta, but this time span 
leaves undetermined the behavior of the low wavenum- 
ber modes over longer time scales. Indeed these modes 
keep growing, as shown in Figure [8j resulting in a steady 
growth of total magnetic energy, shown in Figure [2] 
To understand the long-time dynamics of the system, 
we have performed another simulation, run B, with the 
same physical parameters of run A, but half the orthog- 
onal resolution (Table [1]) , extending the duration up to 
- 11000 r A . 

Figure [9] shows that over longer times the energy of the 
system is prevalently in mode 1, i.e., the largest possible 
scale. But this mode does not grow indefinitely and over 
these much longer time-scales it reaches a statistically 
steady state, fluctuating around its mean value. The 
largest energy fluctuations shown in Figure [9] result in 
energy drops of ~ 2000, that in dimensional units corre- 
spond to a micro-flare with AE ^ 2 x 10 25 erg, releasing 
about twice the amount of energy released by the kink 
instability around t ~ 85 ta in run A (compare with Fig- 
ure [2]). Notice that the kink instability does not appear 
in Figure [9] because the sampling time interval for the 
modes is too long in run B, but it is clearly shown in the 
r.m.s. of the energies (not shown) and in the dissipation 
rate (see inset in Figure HP]). 

Comparing Figure [9] with Figure [TO] where the ohmic 
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dissipation rate J is shown as a function of time, dis- 
plays another interesting result. The large and sharp 
energy drops shown in Figure [9] correspond to large dis- 
sipative peaks in Figure fTQ| e.g., at times t ~ 4750 ta 
and t ~ 10300 ta, but it is also possible to have equally 
large but more gradual energy drops, e.g., between times 
t ~ 6750 ta and t ~ 7750 ta, without a correspond- 
ing single large dissipative peak but rather a cluster of 
smaller peaks. 

In physical space we have already seen in run A that 
initially the inverse cascade corresponds to a perturbed 
magnetic field that occupies an increasingly larger vol- 
ume (Figure 2]) until all the field lines in the box get 
twisted (Figure [5]). In run B we observe that successively, 
once the computational box has been filled with perpen- 
dicular magnetic field, the rising amplitude of modes 1 
corresponds to an increase of the magnetic field inten- 
sity, while the fluctuations in the energy mode are due 
to magnetic reconnect ion events. In fact due to the pe- 
riodic boundary conditions in x and y the same system 
repeats indefinitely along these directions. When the or- 
thogonal magnetic field reaches the boundary it starts 
to interact with the neighboring structures (i.e. with it- 
self coming from the other side). The magnetic energy 
drops in mode 1 correspond to magnetic reconnection 
events that make the system oscillate between the dif- 
ferent possible configurations with energy contained at 
the (large) scales of mode 1 shown in Figure [11] (there is 
no preferred orthogonal direction for the system at this 
scale). 

While the periodic boundary conditions limit the in- 
teractions of large-scale twisted magnetic structures it is 
clearly shown that interaction with such other magnetic 
copies of itself is one of the ways in which the accumu- 
lated energy can be released. Further possibilities and 
the dynamics of these interactions will be the subject of 
future works. 

5. CONCLUSIONS AND DISCUSSION 

In this paper we have investigated the dynamics 
of a closed coronal region driven at its boundary by 
a localized photospheric vortex. Such small vorti- 
cal motions with scales typical of photospheric con- 
vection (~ 1000 km) have been recently observed in 
the p hotosphere ([Brandt et al.l 119881 : iBonet et al J I2QQ8L 
2010J), and can induce releva nt dy namics in the so- 
lar corona (Velli & Liewer 1999; Wede meyer-Bohm et al.l 
120121 : iPanasenco. Martin, fc Vellill2013D . 

A "straightened out" closed region of the solar corona 
is modeled as an elongated Cartesian box where the top 
and bottom plates mimic the photosphere, and the dy- 
namics are integrated with the Reduced M HD equations 
(jKadomtsev fc Pogutsel 11974 iStraussI [l976h . well suited 
for a plasma threaded by a strong axial magnetic field. 

The initial condition consists simply of a uniform ax- 
ial magnetic field. Its field lines are originally straight 
and its footpoints are line-tied at both ends in the top 
and bottom photospheric plates. The photospheric vor- 
tex drags the field lines' footpoints twisting the mag- 
netic field lines (Figure t = 60.57 ta). Even though 
the vortex that we employ is not perfectly circular (Fig- 
ure [TJ Equations (O-©) in the linear stage the field 
lines' tension straightens out in a round shape the or- 
thogonal magnetic field lines (Figure |4j t = 0.61 ta and 



t = 80.64 ta), in this way the Lorentz force vanishes in 
the planes and the system is able to accumulate energy. 
The small departure from a round shape (at its edge) 
of the boundary vortex introduces a small perturbation 
in the coronal field. The system is then unstable to the 
internal kink mode (Figure |4j t = 80.64 ta), and releases 
about 90% of the accumulated energy in a dissipative 
event (Figures [2] and [3]). The energy released in this event 
is of the order of a micro-flare with AE ~ 10 25 erg. 

These resu lts are in agreement with those of 

iLionello et al.l (|1998[ ). that consider similar initial con- 
ditions, performs a refined linear analysis, but does not 
employ a boundary forcing, i.e., the field lines are line- 
tied to a motionless photosphere. Therefore the ini- 
tial linear stage and the development of the kink in- 
stability are in agreement with previous works that 
have always employed large-scale smooth fields with no 
broad-band fluctuations as initial conditions, both in the 
case of field lines line-tied to a motionless photosphere 
(iBaty fc Hevv aerts 1996; Velli ^talMMA ILionello eT~a71 
11993 iBrowning et al.l 2008; H ood et al.ll2QQ9h and with a 
boundary driver (Mik ic et al.lll990l : iGerrard et al.l l2002). 

On the other hand in the solar corona perturbations 
are continuously injected from the lower atmospheric lay- 
ers. Numerical simulations (e.g., Rappaz zo et al.l [2008) 
confirm that especially in closed regions, where waves 
cannot escape toward the interplanetary medium, broad- 
band magnetic fluctuations of the order of a few percent 
of the strong axial magnetic field (not infinitesimal per- 
turbations as classically used in instabilities studies) are 
naturally present. 

Therefore in order to gain a first insight of the coronal 
dynamics when the magnetic field is already structured, 
i.e., there are finite magnetic fluctuations (small but not 
infinitesimal) with small scales and current sheets, we 
continue the simulation after kink instability develops. 
In fact right after kink instability the magnetic energy 
is small (Figure EJ), with b±/Bo ~ 5%, but it is already 
structured with current sheets (Figure HI t = 100.78 ta) 
and a broad band spectrum (Figure [7]). 

The boundary vortex continues to twist the magnetic 
field lines, but in a disordered way (Figure[5j t = 202.15 - 
1211.78 r^). The presence of an already structured mag- 
netic field allows nonlinear dynamics to develop: once 
current sheets and small scales are present, an energy 
cascade continues to feed them, as shown by the energy 
spectra in Figure Therefore current sheets do not dis- 
appear, and the continuous transfer of energy from the 
large to the small scales prevents the field lines to increase 
their twist beyond ~ 180°. The twist remains approxi- 
mately constant in the nonlinear stage as shown in Fig- 
ure [5] (t = 202.15 - 1211.78 ta). Furthermore because 
the current is now concentrated in thin current sheets 
(Figures 3] and [6]) kink instabilities do not develop. 

We had already observed a similar behavior in our pre- 
vious simulations that employed space-filling boundary 
drivers. In particular when the field lines were sheared 
by a ID boundary forcing (Rappazz o. Velli fc Einaudil 
2010) the coronal field was sheared only in the linear 
stage, but after that a multiple tearing instability devel- 
oped and in the coronal field magnetic fluctuations and 
current sheets were formed, the continuous shearing mo- 
tions at the boundary were not able to recreate a sheared 
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Figure 11. Run B: Axial component of the current j (in color) and field lines of the orthogonal magnetic field in the midplane (z = 5) 
at times t ~ 4484 ta and t ~ 4952 ta, i-e., just before and after the dissipative events at time t ~ 4750 ta (Figures [9] and llOJ) . Once the 
orthogonal magnetic field fills the computational box the system oscillates between the different configurations with most of the magnetic 
energy at the large scales, through episodes of magnetic reconnection. 



coronal field and further instabilities were not observed. 

But in the simulations presented in this paper a new 
phenomenon occurs. Although the field lines' twist is ap- 
proximately constant in the nonlinear stage, the volume 
where field lines are twisted increases, and the magnetic 
field acquires larger scales (Figures |4] and [5]). Besides a 
direct cascade that transfers energy from the large to the 
small scales where it is dissipated in current sheets, an 
inverse cascade takes place, transferring energy from the 
injection scale toward larger scales, where no dissipation 
takes place and energy can accumulate. The analysis of 
the magnetic energy modes (Figures [8] and [9]) shows in- 
deed that on long time-scales most of the energy is stored 
at the largest possible scale (mode 1). The inverse cas- 
cade is able to store a significant amount of magnetic 
energy. 

Although magnetic helicity is not defined in RMHD, 
the integral of the magnetic square potential ip is ap- 
proximately conserved (see discussion in last paragraph 
of §4.ip . The inverse cascade of magnetic energy corre- 
sponds also to an inverse cascade of the square potential, 
as clearly shown in Figure |4] where the field lines are the 
contour of i/j (and ip > 0). In future compressible simula- 
tions we expect to observe for magnetic helicity (well de- 
fined in 3D MHD) dynamics similar to those shown here 
for magnetic energy, i.e., an increase of magnetic helicity 
(injected from the boundary) and its inverse cascade, in 
analogy to the the inve rse helicity cascade observed by 
lAmari fc Lucianil (|2000D . 

Because of the periodic boundary conditions along x 
and y the system is virtually repeated along these di- 
rections. When the field lines get twisted in the entire 
computational box, this twisted structure interacts with 
these neighboring twisted structures. This interaction is 
the only condition that limits the growth of magnetic 
energy, giving rise to impulsive magnetic reconnection 
events, that now is not inhibited by the circular topology 
of the orthogonal magnetic field lines of a single struc- 
ture. These events make the system oscillate between the 
many possible configurations with energy in mode 1 (two 



of these are shown in Figure [TT]) . The associated energy 
drops shown in Figure [9] are also in the micro- flare range 
with AE ^ 2 x 10 25 erg, twice the value of the energy 
released initially by the kink instability. 

Although in the presented simulations the generated 
magnetic structures interact only with similar structures 
repeated by the periodic boundary conditions along x 
and y, we can infer that the interaction of a single twisted 
magnetic structure with other magnetic structures can 
give rise to similar release of energy. A more general 
investigations of the interaction between twisted mag- 
netic structures is under way to understand under which 
conditions the interaction leads to energy storage and/or 
release, and to determine quantitatively these properties. 

Previous simulations that employed a space- 
filling; photospheric forcing dRappazzo et all [2008; 
iRappazzo. Velli fc Einaudill2010l ) were not able to accu- 
mulate a significant amount of energy to be successively 
released in micro or larger flares. Those photospheric 
motions, that mimic a uniform and homogeneous 
convection, give instead rise to a basal background 
coronal heating rate in the lower range of the obser- 
vational constr aint ( 1 6 erg cm" 2 s ~ 1 ) and a million 
degree corona (jDahlburg et al.ll2012f ). In the case of a 
space-filling boundary driver we had al so observed that 
the i nverse cascade is inhibited (see Rapp azzo et al.l 
2008, §5.4) for typically strong DC magnetic fields. An 
inverse cascade is pos sible only for weak guide fields (see 
IRappazzo et al.l 12008. §5.4), a condition applicable only 
to limited regions of the corona. 

We conclude that in presence of line- tying and a strong 
guide field, inverse cascade can be a good mechanism to 
store energy, but only if the boundary motion is localized 
in space as the vortex used here, and not space-filling. 
Subsequently the interaction of this magnetic structures 
with others can release the accumulated energy. 

In general photospheric motions will be a superposi- 
tion of approximately homogenous space-filling connec- 
tive motions a nd localized vortical a nd also shearing mo- 
tions (e.g., see Dahlb urg et al.ll2009l . for a localized shear 
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case) . While the space-filling moti ons give rise to a basa l 
background coronal heating (e.g. jRappazzo e~ al. 2008), 
localized motions can give rise to higher impulsive re- 
leases of energy in the micro-flare range and above, con- 
tributing to coronal heating while increasing the tem- 
poral intermittency of the energy deposition and of its 
associated radiative emissions. In future works we will 
consider cases with localized motions superimposed to a 
homogenous space-filling convection-mimicking velocity 
field to determine, among other things, how stronger the 
localized velocity has to be respect to the background 
motions in order to develop dynamics similar to those 
presented in this paper. 

As mentioned in the introduction, highly (and orderly) 
twisted magnetic structures, such &s flux ropes, a r e used 
to initiate solar eruptions (e.g., see iTorok et al.l (120111) 
for a recent application, and the reviews bvlLo^ (|200lh : 
iChenl (|2011f ) for further examples of this model). Kink- 
like instabilities developing in these flux-ropes give rise 
to an explosive dynamics leading to the formation of a 
CME. We have shown that kink-unstable flux ropes are 
not formed in the corona by boundary vortical motions, 
unless a very strong vortex is applied and the coronal 
magnetic fluctuations can then be neglected. Therefore, 
although flux ropes can be formed in the complex dy- 
nami cs in and around a prominence region (jAmari et al.l 
fl999h . given the ubiquitous presence of magnetic fluctu- 
ations in the solar corona, the development of kink-like 
instabilities may be strongly limited. While the dynam- 
ics of the induced CME can be a good approximation, we 
conclude that such models offer a poor model of the ini- 
tiatio n process for whi ch more realistic models are called 
for (|Amari et al.ll2011[ ). 

Generally speaking, in a realistic 3D geometry one 
might expect that the growth of energy in the transverse 
field leads to an inflation and rise of a magnetic loop due 
to the curvature, whi ch we have neglecte d here. This 
effect was included bv lAmari et al.l {1996), who showed 
that twisting the footpoints of a curved flux rope leads 
to its gradual expansion and the system rises to larger 
solar radii. In our simulations the twist does not in- 
crease (the overall field lines twist is limited to 180°), 
remaining roughly constant in the nonlinear stage. It is 
left to future work to understand under which conditions 
such a system, including curvatur e, has dynamics simi- 
lar to those of lAmari et al.l (jl996l), or whether diff erent 
dynamics are possible (see also iGerrard et al1l2004t h and 
how the dynamics develop in a 3D geometry considering 
small or large-scales photospheric vortices. 
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